R语言绘制曼哈顿图示例
其实就是在外观上，类似于这种高楼大厦的风格，换成其它的建筑群，也是很相似的。幸好它已经有名字了，不然我说不定就会直接叫它“shanghaitan plot”，咩哈哈哈
好了不开玩笑了，接正题。
曼哈顿图最常见于全基因组关联分析（Genome-wide association study，GWAS）中，我们常使用它展示基因组中与某种表型显著相关的单核苷酸多态性（SNP）位点或基因型信息，如下来自某项外显子研究，展示了（A）皮肤黑色素瘤和（B）眼黑色素瘤的所有罕见突变位点和肿瘤性状的关联。图中横轴为人22条常染色体，纵轴表示了SNP位点与性状关联的显著性p值，实线表示全基因组统计学显著性阈值，阴影区域表示接近全基因组统计学显著性。
（文献来源：Rare Variant, Gene-Based Association Study of Hereditary Melanoma Using Whole-Exome Sequencing）
非变异位点的关联研究，也可使用曼哈顿图可视化。
例如某项增强子RNA（eRNA）的相关研究，曼哈顿图展示了转录因子（TF）与eRNA关系的分析。图中横轴为人22条常染色体+XY染色体，纵轴表示与每个TF相关的eRNA的百分比，其中红点表示与≥25％的个体eRNA显著相关的假定的主调节转录因子。
（文献来源：Transcriptional landscape and clinical utility of enhancer RNAs for eRNA-targeted therapy in cancer）
或者非人基因组/转录组领域，也仍然适用，毕竟曼哈顿图本身仅仅是一种统计图形而已。比方说这篇微生物组研究中，使用曼哈顿图表示了突变体植株和野生型植株根系微环境中发生显著富集的OTUs（近似理解为微生物物种）。横轴表示不同目（物种分类单元，界门纲目科属种的“目”），纵坐标表示了OTUs差异分析的显著水平，实心点表示该OTU显著富集，点的大小反映了OTUs丰度信息。
（文献来源：Root nodule symbiosis in Lotus japonicus drives the establishment of distinctive rhizosphere, root, and nodule bacterial communities）
在R中，为GWAS分析或绘制曼哈顿图量身打造的包还是非常多的，如qqman等。或者，考虑到曼哈顿图其实就是某种类型的散点图而已，提供作图数据，直接从头绘制也很容易，如使用ggplot2绘制。
本篇分别展示qqman和ggplot2绘制曼哈顿图。将展示两种样式的，如上所示的基因组SNP GWAS的曼哈顿图，以及微生物群落OTUs富集分析的曼哈顿图。
示例作图数据、R脚本等，已上传至百度盘：
https://pan.baidu.com/s/1FGcZ_OYQBGZ6zs4C_AXwzA
GWAS曼哈顿图
qqman包作图
通过qqman包中的作图函数manhattan()，可轻松绘制曼哈顿图。
##qqman 包 manhattan()绘制 GWAS 曼哈顿图library(qqman)
#使用自带的数据集 gwasResults 作示例
#该数据集为人基因组 22 条常染色体中 SNP 位点与某性状的关联程度
data(gwasResults)
head(gwasResults)
#默认作图函数，详情 ?manhattan
#定义 p < 1e-05 为临界显著性，p < 5e-08 为高可信显著性
pdf('GWAS_qq.pdf', width = 12, height = 6)manhattan(gwasResults, col = c('#ff9c00', '#6b6aff'), suggestiveline = -log10(1e-05), genomewideline = -log10(5e-08), annotatePval = 5e-08, annotateTop = FALSE)
dev.off()
ggplot2包作图
我们更换为ggplot2，同样使用上述数据集，绘制曼哈顿图。一个示例如下。
##ggplot2 绘制 GWAS 曼哈顿图library(ggplot2)
#计算染色体刻度坐标
gwasResults$SNP1 <- seq(1, nrow(gwasResults), 1)
gwasResults$CHR <- factor(gwasResults$CHR, levels = unique(gwasResults$CHR))
chr <- aggregate(gwasResults$SNP1, by = list(gwasResults$CHR), FUN = median)
#ggplot2 作图
#定义 p < 1e-05 为临界显著性，p < 5e-08 为高可信显著性
p <- ggplot(gwasResults, aes(SNP1, -log(P, 10))) +
annotate('rect', xmin = 0, xmax = max(gwasResults$SNP1), ymin = -log10(1e-05), ymax = -log10(5e-08), fill = 'gray98') +
geom_hline(yintercept = c(-log10(1e-05), -log10(5e-08)), color = c('blue', 'red'), size = 0.35) +
geom_point(aes(color = CHR), show.legend = FALSE) +
scale_color_manual(values = rep(c('#ff9c00', '#6b6aff'), 11)) +
scale_x_continuous(breaks = chr$x, labels = chr$Group.1, expand = c(0, 0)) +
scale_y_continuous(breaks = seq(1, 9, 2), labels = as.character(seq(1, 9, 2)), expand = c(0, 0), limits = c(0, 9)) +
theme(panel.grid = element_blank(), axis.line = element_line(colour = 'black'), panel.background = element_rect(fill = 'transparent')) +
labs(x = 'Chromosome', y = expression(''~-log[10]~'(P)'))
p
library(ggrepel)
#根据 p 值筛选，将展示子集里的 SNP 标签
gwasResults1 <- subset(gwasResults, P < 5e-08)
#添加 SNP 标签
p1 <- p + geom_text_repel(data = gwasResults1, aes(label = SNP), size = 3,box.padding = unit(0.5, 'lines'), segment.color = 'red', show.legend = FALSE, color = 'red')
p1
ggsave('GWAS_gg.pdf', p1, width = 10, height = 4.5)
#ggsave('GWAS_gg.png', p1, width = 10, height = 4.5)
OTUs富集曼哈顿图
ggplot2包作图
网盘附件“otu_sign.txt”，同样为某16S高通量测序的微生物群落研究，比较了某植株根际环境相对于土壤环境中哪些OTUs发生了显著富集。文件中记录了各OTUs的id、所属的门分类（界门纲目科属种的“门”）、平均相对丰度、差异分析显著性p值以及这些OTUs在根际是否为显著富集的。
然后同样通过曼哈顿图展示这种富集分析结果，就作个和开篇时展示的文献中的那个相似的图吧。
##ggplot2 绘制 OTUs 差异富集曼哈顿图#读取数据
otu_stat <- read.delim('otu_sign.txt', sep = '\t')
#门水平排序，这里直接按首字母排序了
otu_stat <- otu_stat[order(otu_stat$phylum), ]
otu_stat$otu_sort <- 1:nrow(otu_stat)
#其它自定义排序，例如想根据 OTU 数量降序排序
#phylum_num <- phylum_num[order(phylum_num, decreasing = TRUE)]
#otu_stat$phylum <- factor(otu_stat$phylum, levels = names(phylum_num))
#otu_stat <- otu_stat[order(otu_stat$phylum), ]
#otu_stat$otu_sort <- 1:nrow(otu_stat)
#计算 x 轴标签、矩形区块对应的 x 轴位置
phylum_num <- summary(otu_stat$phylum)
phylum_range <- c(0, phylum_num[1])
phylum_name <- phylum_num[1] / 2
for (i in 2:length(phylum_num)) {
phylum_range[i+1] <- phylum_range[i] + phylum_num[i]
phylum_name[i] <- phylum_range[i] + phylum_num[i] / 2
}
#推荐先调整好背景（矩形区块），再绘制散点，这样散点位于矩形区块图层的上面，更清晰
#整体布局，坐标轴、背景线、标签等
p <- ggplot(otu_stat, aes(otu_sort, -log(p_value, 10))) +
labs(x = NULL, y = expression(''~-log[10]~'(P)'), size = 'relative abundance (%)', shape = 'significantly enriched') +
theme(panel.grid = element_blank(), axis.line = element_line(colour = 'black'), panel.background = element_rect(fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) +
scale_x_continuous(breaks = phylum_name, labels = names(phylum_num), expand = c(0, 0)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
#矩形背景
for (i in 1:(length(phylum_range) - 1)) p <- p + annotate('rect', xmin = phylum_range[i], xmax = phylum_range[i+1], ymin = -Inf, ymax = Inf, fill = ifelse(i %% 2 == 0, 'gray95', 'gray85'))
#最后绘制散点：颜色代表门分类，直接使用 ggplot2 默认颜色；大小表示丰度；实心表示富集 OTUs
p <- p +
geom_point(aes(size = abundance, color = phylum, shape = enrich)) +
scale_size(range = c(1, 5))+
scale_shape_manual(limits = c('sign', 'no-sign'), values = c(16, 1)) +
guides(color = 'none') +
geom_hline(yintercept = -log10(0.01), color = 'gray', linetype = 2, size = 1)
p
#输出图片至本地
ggsave('otu.pdf', p, width = 10, height = 5)
#ggsave('otu.png', p, width = 10, height = 5)
后期可以再拿AI、PS等工具修一下。